library(cregg)
library(ggplot2)
library(ggpubr)
rm(list=ls())

#### Load data 

load("data.RData")

#### Formulas for models

f_add <- vote_choice ~ age + prof+career + gender + quality_all
f_add2 <- better_pol ~ age + prof+career + gender + quality_all
f_add3 <- prob_vote ~ age + prof+career + gender + quality_all


#### Compute marginal means and AMCEs

mm_add.vote_choice <- mm(data[data$country!= "Germany" & data$country!="Switzerland",], f_add, id = ~id, family = binomial(link = "logit"))
amces_add.vote_choice <- cj(data[data$country!= "Germany" & data$country!="Switzerland",], f_add, id = ~id)

mm_add.better_pol <- mm(data[data$country!= "Germany" & data$country!="Switzerland",], f_add2, id = ~id, family = binomial(link = "logit"))
amces_add.better_pol <- cj(data[data$country!= "Germany" & data$country!="Switzerland",], f_add2, id = ~id)

mm_add.prob_vote <- mm(data[data$country!= "Germany" & data$country!="Switzerland",], f_add3, id = ~id, family = binomial(link = "logit"))
amces_add.prob_vote <- cj(data[data$country!= "Germany" & data$country!="Switzerland",], f_add3, id = ~id)

#### Create data for Figure A3


dataPlotMM_add.vote_choice <- as.data.frame(list(outcome = c(rep("vote_choice", 21)),
                                                 type = c(rep("MM", 21)),
                                                 level = c("No positive", as.character(mm_add.vote_choice[15,]$level),
                                                           "One positive", as.character(mm_add.vote_choice[16:19,]$level),
                                                           "Two positive", as.character(mm_add.vote_choice[20:25,]$level),
                                                           "Thee positive", as.character(mm_add.vote_choice[26:29,]$level),
                                                           "Four positive", as.character(mm_add.vote_choice[30,]$level)),
                                                 estimate = c(NA, mm_add.vote_choice[15,]$estimate, NA, mm_add.vote_choice[16:19,]$estimate, NA, mm_add.vote_choice[20:25,]$estimate, NA, mm_add.vote_choice[26:29,]$estimate, NA, mm_add.vote_choice[30,]$estimate),
                                                 lower = c(NA, mm_add.vote_choice[15,]$lower, NA, mm_add.vote_choice[16:19,]$lower, NA, mm_add.vote_choice[20:25,]$lower, NA, mm_add.vote_choice[26:29,]$lower, NA, mm_add.vote_choice[30,]$lower),
                                                 upper = c(NA, mm_add.vote_choice[15,]$upper, NA, mm_add.vote_choice[16:19,]$upper, NA, mm_add.vote_choice[20:25,]$upper, NA, mm_add.vote_choice[26:29,]$upper, NA, mm_add.vote_choice[30,]$upper)))

dataPlotamces_add.vote_choice <- as.data.frame(list(outcome = c(rep("vote_choice", 21)),
                                                    type = c(rep("ACME", 21)),
                                                    level = c("No positive", as.character(amces_add.vote_choice[15,]$level),
                                                              "One positive", as.character(amces_add.vote_choice[16:19,]$level),
                                                              "Two positive", as.character(amces_add.vote_choice[20:25,]$level),
                                                              "Thee positive", as.character(amces_add.vote_choice[26:29,]$level),
                                                              "Four positive", as.character(amces_add.vote_choice[30,]$level)),
                                                    estimate = c(NA, amces_add.vote_choice[15,]$estimate, NA, amces_add.vote_choice[16:19,]$estimate, NA, amces_add.vote_choice[20:25,]$estimate, NA, amces_add.vote_choice[26:29,]$estimate, NA, amces_add.vote_choice[30,]$estimate),
                                                    lower = c(NA, amces_add.vote_choice[15,]$lower, NA, amces_add.vote_choice[16:19,]$lower, NA, amces_add.vote_choice[20:25,]$lower, NA, amces_add.vote_choice[26:29,]$lower, NA, amces_add.vote_choice[30,]$lower),
                                                    upper = c(NA, amces_add.vote_choice[15,]$upper, NA, amces_add.vote_choice[16:19,]$upper, NA, amces_add.vote_choice[20:25,]$upper, NA, amces_add.vote_choice[26:29,]$upper, NA, amces_add.vote_choice[30,]$upper)))



dataPlotMM_add.vote_choice$level <- factor(dataPlotMM_add.vote_choice$level, levels= c("No positive", as.character(mm_add.vote_choice[15,]$level),
                                                                                       "One positive", as.character(mm_add.vote_choice[16:19,]$level),
                                                                                       "Two positive", as.character(mm_add.vote_choice[20:25,]$level),
                                                                                       "Thee positive", as.character(mm_add.vote_choice[26:29,]$level),
                                                                                       "Four positive", as.character(mm_add.vote_choice[30,]$level)))
dataPlotamces_add.vote_choice$level <- factor(dataPlotamces_add.vote_choice$level, levels= c("No positive", as.character(mm_add.vote_choice[15,]$level),
                                                                                             "One positive", as.character(mm_add.vote_choice[16:19,]$level),
                                                                                             "Two positive", as.character(mm_add.vote_choice[20:25,]$level),
                                                                                             "Thee positive", as.character(mm_add.vote_choice[26:29,]$level),
                                                                                             "Four positive", as.character(mm_add.vote_choice[30,]$level)))


#### Create data for Figure A4


dataPlotMM_add.better_pol <- as.data.frame(list(outcome = c(rep("better_pol", 21)),
                                                type = c(rep("MM", 21)),
                                                level = c("No positive", as.character(mm_add.better_pol[15,]$level),
                                                          "One positive", as.character(mm_add.better_pol[16:19,]$level),
                                                          "Two positive", as.character(mm_add.better_pol[20:25,]$level),
                                                          "Thee positive", as.character(mm_add.better_pol[26:29,]$level),
                                                          "Four positive", as.character(mm_add.better_pol[30,]$level)),
                                                estimate = c(NA, mm_add.better_pol[15,]$estimate, NA, mm_add.better_pol[16:19,]$estimate, NA, mm_add.better_pol[20:25,]$estimate, NA, mm_add.better_pol[26:29,]$estimate, NA, mm_add.better_pol[30,]$estimate),
                                                lower = c(NA, mm_add.better_pol[15,]$lower, NA, mm_add.better_pol[16:19,]$lower, NA, mm_add.better_pol[20:25,]$lower, NA, mm_add.better_pol[26:29,]$lower, NA, mm_add.better_pol[30,]$lower),
                                                upper = c(NA, mm_add.better_pol[15,]$upper, NA, mm_add.better_pol[16:19,]$upper, NA, mm_add.better_pol[20:25,]$upper, NA, mm_add.better_pol[26:29,]$upper, NA, mm_add.better_pol[30,]$upper)))

dataPlotamces_add.better_pol <- as.data.frame(list(outcome = c(rep("better_pol", 21)),
                                                   type = c(rep("ACME", 21)),
                                                   level = c("No positive", as.character(amces_add.better_pol[15,]$level),
                                                             "One positive", as.character(amces_add.better_pol[16:19,]$level),
                                                             "Two positive", as.character(amces_add.better_pol[20:25,]$level),
                                                             "Thee positive", as.character(amces_add.better_pol[26:29,]$level),
                                                             "Four positive", as.character(amces_add.better_pol[30,]$level)),
                                                   estimate = c(NA, amces_add.better_pol[15,]$estimate, NA, amces_add.better_pol[16:19,]$estimate, NA, amces_add.better_pol[20:25,]$estimate, NA, amces_add.better_pol[26:29,]$estimate, NA, amces_add.better_pol[30,]$estimate),
                                                   lower = c(NA, amces_add.better_pol[15,]$lower, NA, amces_add.better_pol[16:19,]$lower, NA, amces_add.better_pol[20:25,]$lower, NA, amces_add.better_pol[26:29,]$lower, NA, amces_add.better_pol[30,]$lower),
                                                   upper = c(NA, amces_add.better_pol[15,]$upper, NA, amces_add.better_pol[16:19,]$upper, NA, amces_add.better_pol[20:25,]$upper, NA, amces_add.better_pol[26:29,]$upper, NA, amces_add.better_pol[30,]$upper)))



dataPlotMM_add.better_pol$level <- factor(dataPlotMM_add.better_pol$level, levels= c("No positive", as.character(mm_add.better_pol[15,]$level),
                                                                                     "One positive", as.character(mm_add.better_pol[16:19,]$level),
                                                                                     "Two positive", as.character(mm_add.better_pol[20:25,]$level),
                                                                                     "Thee positive", as.character(mm_add.better_pol[26:29,]$level),
                                                                                     "Four positive", as.character(mm_add.better_pol[30,]$level)))
dataPlotamces_add.better_pol$level <- factor(dataPlotamces_add.better_pol$level, levels= c("No positive", as.character(mm_add.better_pol[15,]$level),
                                                                                           "One positive", as.character(mm_add.better_pol[16:19,]$level),
                                                                                           "Two positive", as.character(mm_add.better_pol[20:25,]$level),
                                                                                           "Thee positive", as.character(mm_add.better_pol[26:29,]$level),
                                                                                           "Four positive", as.character(mm_add.better_pol[30,]$level)))


#### Create data for Figure A5


dataPlotMM_add.prob_vote <- as.data.frame(list(outcome = c(rep("prob_vote", 21)),
                                               type = c(rep("MM", 21)),
                                               level = c("No positive", as.character(mm_add.prob_vote[15,]$level),
                                                         "One positive", as.character(mm_add.prob_vote[16:19,]$level),
                                                         "Two positive", as.character(mm_add.prob_vote[20:25,]$level),
                                                         "Thee positive", as.character(mm_add.prob_vote[26:29,]$level),
                                                         "Four positive", as.character(mm_add.prob_vote[30,]$level)),
                                               estimate = c(NA, mm_add.prob_vote[15,]$estimate, NA, mm_add.prob_vote[16:19,]$estimate, NA, mm_add.prob_vote[20:25,]$estimate, NA, mm_add.prob_vote[26:29,]$estimate, NA, mm_add.prob_vote[30,]$estimate),
                                               lower = c(NA, mm_add.prob_vote[15,]$lower, NA, mm_add.prob_vote[16:19,]$lower, NA, mm_add.prob_vote[20:25,]$lower, NA, mm_add.prob_vote[26:29,]$lower, NA, mm_add.prob_vote[30,]$lower),
                                               upper = c(NA, mm_add.prob_vote[15,]$upper, NA, mm_add.prob_vote[16:19,]$upper, NA, mm_add.prob_vote[20:25,]$upper, NA, mm_add.prob_vote[26:29,]$upper, NA, mm_add.prob_vote[30,]$upper)))

dataPlotamces_add.prob_vote <- as.data.frame(list(outcome = c(rep("prob_vote", 21)),
                                                  type = c(rep("ACME", 21)),
                                                  level = c("No positive", as.character(amces_add.prob_vote[15,]$level),
                                                            "One positive", as.character(amces_add.prob_vote[16:19,]$level),
                                                            "Two positive", as.character(amces_add.prob_vote[20:25,]$level),
                                                            "Thee positive", as.character(amces_add.prob_vote[26:29,]$level),
                                                            "Four positive", as.character(amces_add.prob_vote[30,]$level)),
                                                  estimate = c(NA, amces_add.prob_vote[15,]$estimate, NA, amces_add.prob_vote[16:19,]$estimate, NA, amces_add.prob_vote[20:25,]$estimate, NA, amces_add.prob_vote[26:29,]$estimate, NA, amces_add.prob_vote[30,]$estimate),
                                                  lower = c(NA, amces_add.prob_vote[15,]$lower, NA, amces_add.prob_vote[16:19,]$lower, NA, amces_add.prob_vote[20:25,]$lower, NA, amces_add.prob_vote[26:29,]$lower, NA, amces_add.prob_vote[30,]$lower),
                                                  upper = c(NA, amces_add.prob_vote[15,]$upper, NA, amces_add.prob_vote[16:19,]$upper, NA, amces_add.prob_vote[20:25,]$upper, NA, amces_add.prob_vote[26:29,]$upper, NA, amces_add.prob_vote[30,]$upper)))



dataPlotMM_add.prob_vote$level <- factor(dataPlotMM_add.prob_vote$level, levels= c("No positive", as.character(mm_add.prob_vote[15,]$level),
                                                                                   "One positive", as.character(mm_add.prob_vote[16:19,]$level),
                                                                                   "Two positive", as.character(mm_add.prob_vote[20:25,]$level),
                                                                                   "Thee positive", as.character(mm_add.prob_vote[26:29,]$level),
                                                                                   "Four positive", as.character(mm_add.prob_vote[30,]$level)))
dataPlotamces_add.prob_vote$level <- factor(dataPlotamces_add.prob_vote$level, levels= c("No positive", as.character(mm_add.prob_vote[15,]$level),
                                                                                         "One positive", as.character(mm_add.prob_vote[16:19,]$level),
                                                                                         "Two positive", as.character(mm_add.prob_vote[20:25,]$level),
                                                                                         "Thee positive", as.character(mm_add.prob_vote[26:29,]$level),
                                                                                         "Four positive", as.character(mm_add.prob_vote[30,]$level)))


### Plot Figure A3
png("Figures/Figure A3.png", width=2000, height = 3000, res=300)
ggplot(rbind(dataPlotMM_add.vote_choice, dataPlotamces_add.vote_choice))+
  geom_point(aes(x=estimate, y=level))+
  geom_errorbarh(aes(y=level, xmin=lower, xmax=upper), size=0.5, height=0.5)+
  facet_wrap(~type,
             labeller = labeller(
               type = c(`Marginal Means` = "Marginal Means", `Acmes` = "Average Marginal Component Effect")
             ), scales = 'free_x')+
  ylab("")+
  geom_vline(aes(xintercept = c(rep(0.5, 21), rep(0, 21))), linewidth=.1)+
  theme_minimal()+
  scale_y_discrete(limits = rev)+
  theme(legend.position = "bottom", 
        legend.title = element_blank(), 
        panel.border = element_rect(colour = "black", fill=NA, linewidth=.5),
        axis.text.y = element_text(face = c("plain", "bold", 
                                            rep("plain", 4), "bold",
                                            rep("plain", 6), "bold",
                                            rep("plain", 4), "bold",
                                            rep("plain", 1), "bold")))
dev.off()



### Plot Figure A4
png("Figures/Figure A4.png", width=2000, height = 3000, res=300)
ggplot(rbind(dataPlotMM_add.better_pol, dataPlotamces_add.better_pol))+
  geom_point(aes(x=estimate, y=level))+
  geom_errorbarh(aes(y=level, xmin=lower, xmax=upper), size=0.5, height=0.5)+
  facet_wrap(~type,
             labeller = labeller(
               type = c(`Marginal Means` = "Marginal Means", `Acmes` = "Average Marginal Component Effect")
             ), scales = 'free_x')+
  ylab("")+
  geom_vline(aes(xintercept = c(rep(0.5, 21), rep(0, 21))), linewidth=.1)+
  theme_minimal()+
  scale_y_discrete(limits = rev)+
  theme(legend.position = "bottom", 
        legend.title = element_blank(), 
        panel.border = element_rect(colour = "black", fill=NA, linewidth=.5),
        axis.text.y = element_text(face = c("plain", "bold", 
                                            rep("plain", 4), "bold",
                                            rep("plain", 6), "bold",
                                            rep("plain", 4), "bold",
                                            rep("plain", 1), "bold")))
dev.off()




### Plot Figure A5
png("Figures/Figure A5.png", width=2000, height = 3000, res=300)
ggplot(rbind(dataPlotMM_add.prob_vote, dataPlotamces_add.prob_vote))+
  geom_point(aes(x=estimate, y=level))+
  geom_errorbarh(aes(y=level, xmin=lower, xmax=upper), size=0.5, height=0.5)+
  facet_wrap(~type,
             labeller = labeller(
               type = c(`Marginal Means` = "Marginal Means", `Acmes` = "Average Marginal Component Effect")
             ), scales = 'free_x')+
  ylab("")+
  geom_vline(aes(xintercept = c(rep(0.5, 21), rep(0, 21))), linewidth=.1)+
  theme_minimal()+
  scale_y_discrete(limits = rev)+
  theme(legend.position = "bottom", 
        legend.title = element_blank(), 
        panel.border = element_rect(colour = "black", fill=NA, linewidth=.5),
        axis.text.y = element_text(face = c("plain", "bold", 
                                            rep("plain", 4), "bold",
                                            rep("plain", 6), "bold",
                                            rep("plain", 4), "bold",
                                            rep("plain", 1), "bold")))
dev.off()

